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Supercritical  and  transcritical  real-fluid  mixing 
in  diesel  engine  applications 

By  P.  C.  Ma,  L.  Bravof,  and  M.  Ihme 


A  numerical  framework  for  simulating  the  mixing  of  supercritical  and  transcritical  flu¬ 
ids  with  large  density  ratios  is  presented.  The  accurate  and  robust  modeling  of  real-fluid 
effects  on  the  fuel/air  mixing  process  is  critical  in  characterizing  engine  combustion.  To 
describe  this  mixing  process,  the  double-flux  model  is  extended  to  systems  with  real-fluid 
state  equations.  For  this,  an  effective  specific  heat  ratio  is  introduced  to  eliminate  the 
spurious  pressure  oscillations  that  are  caused  by  the  nonlinearity  of  the  real-fluid  equa¬ 
tion  of  state.  High-order  WENO  spatial  discretization  scheme  is  used  for  the  real-fluid 
flow  calculations  to  account  for  the  large  density  gradients  associated  with  these  flows. 
One-dimensional  single  and  multicomponent  advection  test  cases  and  two-dimensional 
supercritical  nitrogen  jet  mixing  and  transcritical  n-dodecane  injection  test  cases  are 
conducted  to  show  the  capability  of  the  present  study  in  application  to  real-fluid  super¬ 
critical  and  transcritical  flows.  Energy  conservation  errors,  generated  due  to  the  double- 
flux  formulation,  are  acceptable. 


1.  Introduction 

A  simple  thermodynamic  analysis  inside  a  diesel  engine  combustion  chamber  will  show 
that  higher  chamber  pressures  will  result  in  higher  energy  output  as  well  as  engine  ef¬ 
ficiency.  Similar  trends  are  also  found  in  other  applications  such  as  rocket  propulsion 
systems  and  gas  turbine  engines.  This  has  motivated  substantial  research  interests  in 
high-pressure  combustion  system.  When  both  of  the  temperature  and  pressure  exceed 
their  critical  values,  the  thermodynamic  state  is  called  a  supercritical  state.  In  diesel 
engines,  typically  a  subcritical  liquid  fuel  is  injected  into  a  supercritical  ambient  gas  and 
the  fuel  jet  will  be  heated  to  supercritical  temperature  before  combustion.  This  process 
is  often  referred  to  as  a  “transcritical”  injection.  At  elevated  pressure,  the  mixture  prop¬ 
erties  exhibit  liquid-like  densities,  gas-like  diffusivities,  and  surface  tension  and  enthalpy 
of  vaporization  approach  zero  (Yang  2000).  This  phenomenon  has  been  shown  by  recent 
experimental  works  (Mayer  et  al.  2000;  Oschwald  et  al.  2006;  Chehroudi  2012).  Recent 
research  by  Dahrns  &  Oefelein  (2013)  also  provided  conceptual  insights  into  diesel  in¬ 
jection  transition  processes  at  high  pressures.  However,  these  complex  processes  are  still 
not  well  understood  either  experimentally  or  numerically.  Therefore,  to  provide  insight 
into  the  high-pressure  combustion  systems,  accurate  and  robust  simulation  tools  are  re¬ 
quired  for  the  characterization  of  supercritical  and  transcritical  flows  where  large  density 
gradients  and  thermodynamic  anomalies  occur. 

Several  difficulties  require  consideration  for  the  modeling  of  supercritical  fluid  mixing 
numerical  simulations.  One  of  them  is  the  large  density  gradients  caused  by  the  liquid-like 
density  in  the  dense  fluid  region  and  the  gas-like  density  in  the  supercritical  ambient  gas 
region.  A  typical  treatment  in  the  literature  is  to  introduce  artificial  dissipation  locally 
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where  large  density  gradients  occur  (Scllc  &  Schmitt  2010;  Ruiz  2012;  Terashima  & 
Koshi  2012).  Our  previous  work  (Hickey  et  al.  2013;  Hickey  &  Ihme  2013)  implemented 
a  hybrid  central-ENO  scheme  in  analogy  to  shock-capturing  schemes  to  locally  switch  to 
second-order  ENO  scheme  to  stabilize  regions  of  large  density  gradients. 

Another  obstacle  is  spurious  pressure  oscillations  that  are  generated  when  a  fully  con¬ 
servative  scheme  is  adopted,  due  to  the  high  nonlinearity  of  the  real-fluid  equation  of  state 
(RF-EoS),  which  is  similar  to  that  of  multicomponent  compressible  ideal  gas  flow  calcu¬ 
lations  (Abgrall  &  Kami  2001).  To  prevent  this,  Ruiz  (2012)  derived  a  quasi-conservative 
scheme  by  relating  the  artificial  dissipation  terms  in  the  mass,  momentum  and  energy 
conservation  equations  and  setting  the  pressure  differential  to  be  zero.  Terashima  &  Koshi 
(2012)  solved  the  pressure  equation  instead  of  the  energy  equation  based  on  the  work 
of  Abgrall  (1996).  For  ideal  gas,  Johnsen  &  Ham  (2012)  adopted  a  scheme  in  which  an 
auxiliary  advection  equation  for  the  specific  heat  ratio  was  solved.  Although  this  method 
is  suitable  for  well-defined  interfacial  flows  with  inert  species,  it  is  not  applicable  for  ther¬ 
modynamically  complex  supercritical  flows  in  which  the  thermodynamic  properties  are 
dependent  on  temperature  and  species  compositions.  Another  approach  was  proposed  by 
Abgrall  &  Kami  (2001)  and  Billet  &  Abgrall  (2003)  and  later  extended  by  Houim  &  Kuo 
(2011).  This  method  is  referred  to  the  double- flux  model  for  ideal  gas  flows.  However, 
the  double- flux  model  has  so  far  not  been  extended  to  RF-EoS. 

In  this  study,  we  propose  a  numerical  scheme  for  the  simulation  of  real  fluids  mixing  at 
supercritical  and  transcritical  conditions.  To  this  end,  the  double-flux  model  is  extended 
to  the  real-fluid  conditions  by  introducing  the  effective  specific  heat  ratio.  Higlr-oder  ac¬ 
curate  spatial  discretization,  a  WENO  scheme  (Jiang  &  Shu  1996)  and  the  corresponding 
Riemann  solver  were  applied  in  the  real-fluid  context  to  deal  with  large  density  gradients. 
Several  test  cases  in  one  and  two  dimensions  are  conducted  to  demonstrate  the  capability 
of  this  scheme  to  deal  with  complex  thermodynamics  while  maintaining  robustness  and 
high-order  accuracy. 


2.  Mathematical  formulation 


2.1.  Governing  equations  and  equation  of  state 


The  multi-component  non-reacting  Navier-Stokes  equations  in  fully  conservative  form  is 
the  working  set  of  the  current  study.  A  RF-EoS  is  adopted  to  account  for  the  consistent 
thermodynamics  in  supercritical  and  transcritical  regimes.  For  numerical  efficiency,  the 
state  equation  should  be  simple  enough  to  allow  for  a  rapid  computation  of  the  thermo¬ 
dynamic  states  while  simultaneously  accounting  for  the  full  complexity  of  the  nonlinear 
relationship  among  all  state  quantities.  The  cubic  Peng-Robinson  equation  of  state  (Peng 
&  Robinson  1976)  was  chosen  since  it  provides  a  computationally  efficient  and  accurate 
state  representation  near  the  critical  point  (Miller  et  al.  2001).  The  cubic  state  equation 
is  written  as, 


1ZT  Am 

V-Bm  ~  V2  +  2 VBm+Bl  ’ 


(2.1) 


where  V  and  1Z  represent,  respectively,  the  molar  volume  and  universal  gas  constant. 
The  coefficients  Am  and  Brn  are  evaluated  using  the  critical  properties  and  temperature 
of  the  fluid.  As  these  equations  are  developed  for  pure  fluids,  mixing  rules  are  applied 
to  treat  multi-species  mixtures.  A  comparison  of  the  Peng-Robinson  and  Soave-Redlicli- 
Kwong  (Soave  1972)  equations  of  state  against  NIST  database  for  n-dodecane  under  the 
“Spray  A”  conditions  (Pickett  &  Bruneaux  2011)  is  presented  in  Figure  1.  The  method 
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Figure  1 .  Density  of  n-dodecane  at  6  MPa  predicted  by  Peng- Robinson  ( - )  and 

Soave-Redlich-Kwong  (-  -)  equations  of  state  agains  NIST  database  (o)  and  ideal  gas  law 

( - ).  Vertical  line  shows  the  critical  temperature  of  n-dodecane  ( Tc  =  658.2  K).  The  critical 

pressure  of  n-dodecane  pc  =  1.80  MPa. 


by  Chung  et  al.  (1988)  is  used  to  evaluate  the  dynamic  viscosity  of  the  mixture.  Details 
on  the  mixing  rules  and  transport  models  can  be  found  in  Poling  et  al.  (2001). 


2.2.  Double-flux  model 

To  prevent  spurious  pressure  oscillations,  the  double-flux  model  is  extended  to  real-fluid 
simulations.  The  double-flux  model  was  proposed  by  Abgrall  &  Kami  (2001),  for  the 
computations  of  compressible  multicomponent  flow.  Analysis  on  the  Euler  system  shows 
that  when  a  fully  conservative  scheme  is  adopted  for  multicomponent  calorically  perfect 
gas,  the  pressure  equilibrium  is  not  be  preserved  and  spurious  oscillations  occur.  The 
occurrence  of  these  instabilities  are  attributed  to  the  different  values  of  the  specific  heat 
ratio,  7,  since  for  calorically  perfect  gas,  the  coupling  between  total  energy  and  pressure 
is, 


E  =  pe+  -putUi 


P 

7-1 


1 

—  piLiUi  , 


(2.2) 


where  E  and  e  are  total  energy  and  specific  internal  energy,  respectively.  For  multicompo¬ 
nent  cases,  7  is  a  nonlinear  function  of  gas  compositions.  Note  that  for  single  component 
calorically  perfect  systems,  in  which  7  remains  constant,  no  pressure  oscillation  are  ob¬ 
served  due  to  the  linear  property  of  Eq.  (2.2). 

The  simple  single-fluid  method  proposed  by  Abgrall  &  Kami  (2001)  assumes  that  ad¬ 
jacent  cells  have  the  same  7  value  as  that  of  the  cell  itself,  when  computing  the  flux.  Note 
that  by  fixing  7  the  scheme  is  no  longer  fully  conservative.  Energy  conservation  errors  are 
expected.  Abgrall  &  Kami  (2001)  show  that  this  error  is  proportional  to  the  difference  in 
7  and  converges  with  increasing  resolution.  In  combustion  applications  or  supercritical 
jet  mixing  cases,  mach  number  is  typically  not  very  high  so  the  small  conservation  errors 
are  not  expected  to  cause  a  problem  (Lv  &  Ihrne  2014).  The  nonlinearity  between  the 
coupling  of  the  total  energy  and  the  pressure  can  arise  from  the  temperature  dependent 
7  typically  in  combustion  applications.  Billet  &  Abgrall  (2003)  extended  the  original 
double-flux  model  to  account  for  variable  7  values  due  to  variations  in  both  temperature 
and  gas  compositions. 

For  the  RF-EoS,  the  simple  relation  between  specific  internal  energy  and  pressure 
through  7  is  no  longer  valid  and  instead  they  are  coupled  by  an  integral  equation.  For 
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Case  Pressure 


0.25  m  <  x  <  0.75  m 


x  <  0.25  m  or  x  >  0.75  m  Density  ratio 


A 

B 

C 


4  MPa 

5  MPa 

6  MPa 


N2,  440.0  kg/m3,  126.9  K 
Ho,  25.8  kg/m3,  52.8  K 
n-Dodecane,  642.0  kg/m3,  363 


N2,  44.5  kg/m3,  300  K 
N2,  50.8  kg/m3,  332.2  K 
K  N2,  22.1  kg/rn3,  900  K 


10 

2 

29 


Table  1.  Initial  conditions  for  the  one-dimensional  test  cases.  Domain  length  is  a;  £  [0, 1]  m. 


single  component  cases,  the  expression  of  the  specific  internal  energy  in  real  fluids  is, 


j-Ti 

e(T\,  p\)  =  eref  +  /  cv(T)dT  - 1- 

"  ^"ref 


(2.3) 


where  eref  is  the  reference  ideal  gas  specihc  internal  energy  value  at  Tref,  and  is  the 
ideal  gas  specific  heat  capacity  which  is  evaluated  using  the  NASA  polynomials  in  the 
current  study.  The  last  integral  term  in  Eq.  (2.3)  can  be  evaluated  analytically  based  on 
the  cubic  equation  of  state  used.  Note  that  here  the  specific  internal  energy  is  a  function 
of  both  temperature  and  density,  which  is  different  from  that  in  the  ideal  gas  case.  To 
extend  the  double-flux  model  to  the  RF-EoS,  we  propose  a  simple  effective  specific  heat 
ratio,  7,  which  is  defined  as  follows, 

7  =  1  +  —  -  (2.4) 

pe 

This  definition  is  derived  by  relating  the  specific  internal  energy  and  pressure  through 
pe  =  p/{j  —  1)  in  analogy  to  the  formula  for  specific  internal  energy  in  the  ideal  gas 
case  as  in  Eq.  (2.2).  Through  this  effective  specific  heat  ratio,  similar  procedures  of  the 
double-flux  model  can  be  applied  ensuring  that  the  pressure  equilibrium  will  be  preserved 
as  shown  in  the  following  one-dimensional  test  cases. 


2.3.  Spatial  discretization  and  time  advancement 

A  finite  volume  spatial  discretization  scheme  is  adopted  in  the  current  study.  A  WENO 
reconstruction  scheme  was  used  to  evaluate  the  left  and  right  state  values  for  the  Rie- 
mann  solver.  This  provides  high-oder  accuracy  and  essentially  non-oscillatory  behavior. 
Primitive  variables  are  interpolated  onto  the  cell  faces.  Other  choices  of  interpolation 
variables  are  also  possible,  such  as  conservative  variables  and  characteristic  variables, 
and  this  will  be  evaluated  in  future  work.  The  procedure  for  converting  primitive  vari¬ 
ables  to  conservative  variables  and  vice  versa  can  be  found  in  Hickey  et  al.  (2013)  and 
Hickey  &  Ihme  (2013).  A  TVD  RK3  (Gottlieb  &  Shu  1998)  time  advancement  scheme  is 
adopted  in  the  current  study. 

In  the  current  study,  the  HLLC  approximate  Riemann  solver  is  used  for  the  solution 
of  the  Riemann  problem.  Several  approaches  are  possible  to  estimate  the  wave  speeds 
(Batten  et  al.  1997),  and  here  the  minimum  and  maximum  of  the  left  and  right  going 
characteristic  wave  speeds  are  used. 

Similar  to  that  of  Houim  &  Kuo  (2011),  the  procedure  for  time  advancement  from  time 
tn  to  time  tn+1  using  the  double-flux  model  in  current  study,  is  summarized  as  follows, 

1.  Compute  and  store  effective  specific  heat  ratio  'yn  according  to  Eq.  (2.4)  for  each  cell. 

2.  Reconstruct  the  face  values  of  primitive  variables  based  on  the  chosen  spatial  dis¬ 
cretization  scheme  and  construct  the  double- flux,  based  on  the  RF-EoS. 
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Figure  2.  Simulation  results  after  one  period  at  t  =  0.02  s  for  case  A.  Clockwise  density, 

pressure,  velocity  and  temperature  profiles  respectively.  - ,  exact  solution;  — o — ,  WEN05; 

— 0 — ,  WEN03;  — □—  ,  first-order  upwinding. 


3.  Update  conservative  variables  using  the  double-flux  model  and  the  time  marching 
method. 

4.  Update  primitive  variables.  For  pressure,  p  =  (7"  —  1  )(E  —  \putUi). 

5.  Repeat  steps  2-4  for  each  stage  of  the  time  integration  scheme. 

6.  After  all  stages  of  the  time  integration  scheme  are  done,  update  conservative  variables 
from  primitive  variables. 


3.  One-dimensional  test  cases 

Three  one-dimensional  advection  configurations  have  been  selected  to  evaluate  the 
performance  of  the  proposed  algorithm.  Note  that  Euler  equations  are  solved  here, 
enabling  a  direct  comparison  with  the  analytical  solutions.  The  initial  conditions  of 
the  three  test  cases  are  listed  in  Table  1.  The  pressure  for  all  cases  are  all  above 
the  critical  pressure  for  all  species  considered  (pCjn2  =  3.40  MPa,  pCyu2  =  1.28  MPa, 
Pc,n-Dodecane  =  1.80 MPa).  Case  A  is  a  supercritical  case  under  the  same  conditions  as 
the  case  3  of  Mayer  et  al.  (2003).  Cases  B  and  C  are  transcritical  cases,  meaning  that 
the  temperature  of  the  dense  fluid  is  below  its  critical  temperature  (Tc,n2  =  126.2  K, 
TCih2  =  32.9 K,  TCin. Dodecane  =  658.2  K).  Case  C  corresponds  to  “Spray  A”  condition 
(Pickett  &  Bruneaux  2011)  as  in  diesel  engine  applications.  The  computational  domain 
is  a;  £  [0, 1]  m  and  for  all  computations  51  grid  points  are  used,  giving  a  uniform  grid 
spacing  of  Aa;  =  0.02  m.  Periodic  boundary  conditions  are  applied.  For  all  computations, 
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Figure  3.  Simulation  results  after  one  period  at  t  =  0.02  s  for  case  B.  Density  profile  on  the 
left  and  hydrogen  mass  fraction  on  the  right.  Same  notations  as  in  Figure  2. 


Figure  4.  Simulation  results  after  one  period  at  t  =  0.02  s  for  case  C.  Density  profile  on  the 
left  and  n-dodecane  mass  fraction  on  the  right.  Same  notations  as  in  Figure  2. 


the  CFL  number  is  set  to  be  1.0.  The  advection  velocity  is  50m/s  and  one  period  is  com¬ 
puted  until  t  =  0.02  s.  Three  spatial  discretization  scheme  are  tested,  5th-order  WENO, 
3rd-order  WENO  and  first-order  upwinding. 

Figure  2  shows  the  results  for  density,  pressure,  velocity  and  temperature  profiles  for 
case  A.  The  present  numerical  scheme  preserves  the  pressure  and  velocity  equilibrium 
as  expected  and  no  spurious  oscillations  are  generated.  Note  that  for  this  cryogenic 
case,  a  fully  conservative  scheme  will  fail  due  to  occurrence  of  the  spurious  pressure 
oscillation  unless  an  extremely  large  number  of  grid  points  is  used  to  increase  the  spatial 
resolution.  The  density  and  temperature  profiles  do  not  exhibit  any  visible  oscillations, 
demonstrating  the  capability  of  the  WENO  scheme  and  the  HLLC  Riemann  solver  in 
the  context  of  real-fluid  simulations.  The  first-order  upwinding  scheme  is  too  diffusive 
in  this  case.  The  energy  conservation  errors  are  evaluated  for  all  cases,  and  they  remain 
below  1%. 

Figures  3  and  4  show  the  results  for  density  and  mass  fraction  profiles  for  cases  B  and  C. 
These  two  cases  demonstrate  the  capability  of  the  newly  developed  scheme  for  application 
to  multicomponent  real-fluid  flows.  Again,  no  spurious  oscillations  of  pressure  or  velocity 
are  generated  (results  not  shown  here) .  And  the  density  and  mass  fraction  profiles  do  not 
show  any  visible  fluctuations.  The  behavior  of  different  spatial  discretization  schemes  is 
similar  to  that  for  case  A. 
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Figure  5.  Simulation  results  of  the  planar  cryogenic  nitrogen  jet  mixing  case  at  £  =  7.90  x  1CP4  s. 
Density  profile  at  top  and  pressure  profile  at  bottom.  Contour  ranges  not  adjusted.  Results  in 
full  computation  domain  (30 h  x  16h  where  h  =  1  mm)  are  shown. 

4.  Two-dimensional  test  cases 

4.1.  Cryogenic  nitrogen  jet  mixing  case 

A  two-dimensional  cryogenic  planar  jet  mixing  case  is  carried  out  to  demonstrate  the 
multidimensional  capability  of  the  current  numerical  scheme.  A  cryogenic  nitrogen  jet 
(126.9  K)  is  injected  into  a  high  temperature  (300  K)  nitrogen  vessel.  The  ambient 
pressure  is  initially  set  to  be  constant  at  4  MPa.  This  case  is  a  supercritical  mixing 
case  and  similar  to  the  condition  of  the  case  A  (see  Table  1)  and  case  3  of  Mayer  et  al. 
(2003).  Note  that  here  the  Navier-Stokes  equations  are  solved  including  the  viscous  terms. 
The  height  of  the  jet  exit  nozzle  is  1.0  mm  and  the  jet  velocity  is  set  to  be  50  m/s. 
A  domain  of  30 h  x  16 h  is  used,  where  h  is  the  jet  height.  A  uniform  mesh  in  both 
directions  is  used,  which  has  a  minimum  spacing  of  0.04/i.  The  inlet  condition  of  the 
jet  is  a  plug  flow  with  a  top-liat  velocity  profile.  Pressure  outlet  boundary  conditions 
are  applied  at  the  top,  bottom  and  right  boundaries.  Adiabatic  wall  is  applied  to  the 
left  wall.  No  subgrid  scale  model  is  used  in  this  case.  Note  that  the  mesh  resolution 
is  not  sufficient  to  resolve  all  scales  in  this  flow,  since  the  purpose  of  this  test  case  is 
to  demonstrate  the  current  numerical  scheme.  The  solver  used  for  the  two-dimensional 
cases  is  in-house,  unstructured,  finite-volume  code,  CharLESx ,  in  which  the  double-flux 
model  is  implemented.  A  hybrid  central-ENO  spatial  discretization  scheme  is  adopted 
here.  The  density-based  sensor  switches  between  a  second-order  ENO  and  a  fourth-order 
central  scheme  to  suppress  spurious  oscillations  generated  due  to  large  density  gradients. 
Details  can  be  found  in  Hickey  et  al.  (2013)  and  Hickey  &  Ihrne  (2013). 

Figure  5  shows  the  density  and  pressure  fields  of  the  nitrogen  jet  mixing  case.  The 
pressure  field  does  not  show  any  spurious  oscillations  in  this  viscous  case.  Note  that  for 
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Figure  6.  Simulation  results  of  the  transcritical  n-dodecane  injection  case  at  time  60  ps, 
120  ps  and  180  ps  respectively  from  top  to  bottom.  Left  column  shows  density  fields  and 
right  column  shows  n-dodecane  mass  fraction  fields.  Only  portion  of  the  computation  domain 
([0,40]  x  [—8,  8]  mm)  is  shown. 

the  viscous  case,  the  pressure  equilibrium  from  the  Euler  system  is  not  expected.  For 
fully  conservative  schemes  as  used  by  Hickey  et  al.  (2013)  and  Hickey  &  Ihme  (2013), 
very  fine  mesh  and  more  diffusive  schemes  are  needed  to  suppress  spurious  oscillations. 
In  contrast,  the  present  scheme  provides  stable  solutions  for  pressure  despite  the  coarse 
mesh. 


4.2.  Transcritical  n-dodecane  injection  case 

With  the  success  on  the  cryogenic  nitrogen  mixing  case,  the  present  numerical  scheme 
is  then  tested  at  conditions  that  are  more  representative  for  diesel  engine  applications. 
For  this,  the  injection  of  n-dodecane  at  500  K  into  a  vessel  of  nitrogen  at  900  K  is 
considered.  The  pressure  is  initially  set  to  be  a  constant  at  6  Mpa.  The  condition  is 
similar  to  the  “Spray  A”  condition  (Pickett  &  Bruneaux  2011).  Adiabatic  wall  boundary 
conditions  are  applied  to  the  left,  top  and  bottom  walls.  The  pressure  outlet  is  applied 
at  the  right  outlet.  Plug  flow  is  applied  at  the  inlet.  The  injection  velocity  is  set  to  be 
650  m/s  based  on  the  mass  flow  rate  measurements.  Navier-Stokes  equations  are  solved 
in  the  CharLESx  solver  with  the  hybrid  central-ENO  scheme  and  the  two-dimension 
computational  domain  of  800/i  x  400/i  (with  h  =  0.09  mm  being  the  injector  height) 
is  used  to  eliminate  boundary  effects.  The  minimum  mesh  spacing  is  4.5  pm  and  total 
number  of  grid  points  is  about  0.2  million.  No  subgrid  scale  model  is  applied,  since  the 
purpose  of  this  test  case  is  not  to  resolve  all  scales  in  the  flow  but  to  demonstrate  the 
capability  of  the  current  scheme  to  diesel  engine  applications. 

Figure  6  shows  simulation  results  for  density  and  n-dodecane  mass  fraction.  During  the 
simulation,  no  spurious  pressure  oscillations  were  generated.  Substantially  high  resolution 
is  needed  for  a  fully  conservative  scheme  to  simulate  similar  conditions  without  significant 
pressure  oscillations.  Temperature  spikes  were  observed  in  the  simulation  results  (not 
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shown  here),  and  this  may  be  solved  by  interpolating  the  characteristic  variable  for  the 
reconstruction  (Houim  &  Kuo  2011).  This  will  be  addressed  as  a  part  of  the  future  work. 
We  have  demonstrated  the  overall  performance  of  the  herein  developed  numerical  scheme 
by  considering  one-  and  two-dimensional  test  cases  involving  density  ratio  in  excess  of 
25  and  the  extension  to  three-dimensional  configurations  is  the  focus  of  future  work. 


5.  Conclusions 

A  numerical  framework  for  simulating  supercritical  and  transcritical  fluid  mixing  with 
large  density  ratios  is  presented.  The  double-flux  model  is  extended  to  RF-EoS  by  in¬ 
troducing  an  effective  specific  heat  ratio  to  eliminate  spurious  pressure  oscillations  that 
are  caused  by  the  highly  nonlinear  RF-EoS.  A  high-order  WENO  spatial  discretiza¬ 
tion  scheme  is  used  along  with  a  TVD  RK3  time  integration  scheme  to  account  for 
large  density  gradients  associated  with  these  flows.  One-dimensional  advection  test  cases 
demonstrate  the  superior  performance  of  the  present  numerical  scheme  over  fully  con¬ 
servative  schemes.  The  conservation  errors  generated  due  to  the  double-flux  formulation 
are  acceptable.  No  spurious  pressure  oscillations  are  observed  for  all  three  test  cases  and 
no  oscillations  in  density  field  are  observed.  Two-dimensional  test  cases  of  supercritical 
nitrogen  jet  mixing  and  the  transcritical  n-dodecane  injection  show  the  capability  of  the 
present  study  to  multidimensional  flows.  No  spurious  pressure  oscillations  are  observed 
even  on  relatively  coarse  meshes.  Future  work  includes  the  application  of  the  method  to 
three-dimensional  configurations  and  the  validation  against  experimental  measurements. 
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